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A major challenge in realizing antiferromagnetic (AF) and superfluid phases in optical lattices 
is the ability to cool fermions. We determine the equation of state for the 3D repulsive Fermi- 
Hubbard model as a function of the chemical potential, temperature and repulsion using unbiased 
determinant al quantum Monte Carlo methods, and we then use the local density approximation to 
model a harmonic trap. We show that increasing repulsion leads to cooling, but only in a trap, due 
to the redistribution of entropy from the center to the metallic wings. Thus, even when the average 
entropy per particle is larger than that required for antiferromagnetism in the homogeneous system, 
the trap enables the formation of an AF Mott phase. 

PACS numbers: 71.10.Fd, 37.10.Jk, 71.27. +a 



Introduction: One of the most exciting themes in 
condensed matter physics is how complex states of mat- 
ter emerge from simple Hamiltonians. In particular, the 
repulsive Fermi-Hubbard model gives rise to a rich vari- 
ety of behavior, including a Mott insulating regime, an 
antiferromagnetically ordered Neel state, and possibly a 
"high-temperature" d-wave superfluid. 

Cold atomic gases are unique in being clean and tun- 
able systems that offer tremendous promise for explor- 
ing such Hamiltonians. The Fermi-Hubbard model can 
be emulated using an optical lattice with two hyper- 
fine species of fermions [2]. Several experimental feats 
have already been accomplished: the observation of sharp 
Fermi surfaces for free fermions in an optical lattice [3], 
and of the Mott insulating regime for repulsively inter- 
acting fermions [4, 5]. The next step in this quest is to 
go to even lower temperatures, where the local moments 
order to form a Neel antiferromagnet. 

In this Letter we present an adiabatic cooling protocol 
for trapped systems, which we expect to play an impor- 
tant role in the race for finding antiferromagnetism in 
the repulsive Hubbard model and for opening the door 
toward the search for the d-wave superfluid state. We 
first calculate the thermodynamics of a homogeneous sys- 
tem using unbiased determinantal quantum Monte Carlo 
(DQMC) as a function of filling and temperature, access- 
ing both paramagnetic and AF phases. At half-filling, 
this allows us to obtain the entropy down to T = O.lt 
(see Fig. 1(b)), well below the maximum Neel tempera- 
ture Tn ~ 0.36t [6], and also well below the temperatures 
accessed by recent cluster studies [1]. 

We next use the local density approximation to treat 
the effect of a harmonic trap. We demonstrate that in- 
creasing the repulsion U adiabatically leads to substan- 
tial cooling, but only in the presence of the trap (see 
Fig. 2). During this process, the cloud expands and en- 
tropy gets redistributed from the center to the metal- 
lic wings. Even though the average entropy per particle 
S/N « is higher than the critical entropy of the 
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FIG. 1: (a) Energy per site of homogeneous system 

at half- filling and U/t = 8, calculated using DQMC down to 
T/t = 0.1. Statistical error bars are smaller than symbols. 
The solid curve is the entropy extrapolated to L = oo and 
dr = (details in supplement), (b) Entropy per site ob- 
tained by integrating down from T = oo, showing a shoulder 
at the Mott scale Tyiott — U and a distinct feature at the Neel 
temperature Tn ~ 0.36£ due to critical fluctuations. Errors 
in E/t and s/ks are both about 0.02. DCA results [1] are 
shown for comparison. 
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homogeneous system (OAks at U/t = 8), it is neverthe- 
less possible to generate an AF state at the center (see 
Fig. 3). 

Model and methods: We consider the Fermi- 
Hubbard Hamiltonian, 

H = - t ( c i*c r >* + c Lc r >a) + u ^2n rt n ri 

(rr')a r 

+ £(y t r 2 -^)(n r -l), (1) 

r 

in which r labels a site (or well) of a 3D cubic optical 
lattice, a =t or ^ corresponds to two hyperfine states, 
t is the nearest-neighbor hopping amplitude, U is the 
on-site interaction energy, c rcr is the fermion destruc- 
tion operator at site r with spin <j, and n rcr = cl a c ra 
with n r = ^2 a n ra . The curvature V t = \muj\d? de- 
scribes harmonic confinement with trap frequency ojo/2tt, 
fermion mass m, and lattice spacing d. The chemical po- 
tential ji controls the average density. The parameters 
t and U can be directly related [8] to the lattice depth, 
set by the laser intensity, and to the interatomic interac- 
tion tuned by a Feshbach resonance. This Hamiltonian 
is valid in the regime where only a single band is pop- 
ulated in the optical lattice. Following Ref. 5 we define 
the characteristic trap energy E t = ^(3A^/87r) 2//3 . 

We calculate the density p, energy density E, double 
occupancy D = (n r ^n r |), and spin correlations for a ho- 
mogeneous system (V t = 0) as a function of /i, T, and 
U/t using DQMC simulations [9, 10]. 

Half-filling: We first focus on the homogeneous case 
at half-filling (p = 0) and U/t = 8, where the Neel tem- 
perature T N /t = 0.36 is highest [6]. At p = DQMC 
is free of the fermion sign problem and we can access 
low temperatures down to T = O.lt, well into the AF 
phase. We perform extrapolation on E(T) to the limit 
of zero imaginary-time discretization (Sr = 0) and infi- 
nite system size (L 3 = oo), as described in detail in the 
supplement. The high statistical accuracy of the DQMC 
data even reveals critical fluctuations near Tn- 

We obtain the ground state energy Eo/t = —0.74(2) 
and the correct low-temperature behavior (E ~ T 4 ) ex- 
pected for an antiferromagnet with linearly dispersing 
spin waves. The results are shown in Fig. 1(a). Inte- 
grating E(T) down from infinite temperature, we deter- 
mine the entropy per site using s(T) = In 4 + E /T — 
J^dTE/T 2 . Our results agree with extrapolated results 
from the dynamical cluster approximation (DC A) [1], 
available only in the paramagnetic phase. 

We see from Fig. 1(b) that as the temperature is re- 
duced below U — 8t, the entropy per site s/ks decreases 
from ln(4) to ln(2), due to suppression of double occu- 
pancy below the Mott scale for charge fluctuations. At 
Tn the critical entropy is sn/^b ~ 0.4/c^. Our DQMC 
results show a steep drop in entropy below Tn resulting 
from spin ordering. 



In Fig. 2(a), we show constant-entropy curves in the 
(T, U) plane at half-filling. We also plot the Neel tem- 
perature as a function of U obtained from previous 
QMC simulations [6] together with its asymptotic forms 
at weak and strong coupling. The dashed curve is 
0.282Tmf (U/t) where the mean-field result is given by 
2/U = tanh(2e/ e /TMF)/e/e and the suppression factor 
0.282 arises from 0((U/t) 2 ) vertex corrections [11, 12]. 
The dotted curve shows the strong-coupling Heisenberg 
limit result 3.78t 2 /U [13]. 

Away from half- filling: We next compute the equa- 
tion of state p(p) of the homogeneous system away from 
half-filling, as this will be needed to study the effect 
of a trap. We now obtain the entropy by integrat- 
ing along an isotherm from the empty lattice, s(p) = 
f-oed/J* (ds/dfjL)T, making use of the Maxwell relation 
(ds/dp)r — (dp/dT)^ where (dp/dT)^ is evaluated us- 
ing a finite difference scheme. This gives results (indi- 
cated by symbols labelled "J d/i" in Fig. 1(b)) consistent 
with integration of E(T) as described above. 

We model the trap using the local density approxi- 
mation (LDA), in which local observables are given by 
their homogeneous values evaluated at a chemical poten- 
tial p(r) = /io — V t r 2 . The chemical potential at the trap 
center /iq is determined from the total fermion number 
TV = J °° dr Airr 2 p(p(r)) . We obtain density, entropy, 
and local spin correlation profiles such as those in Figs. 3 
and 4, from which we can deduce a route to achieving 
cooling in optical lattices. 

Cooling: Note the contrast between the constant- 
entropy curves in the homogeneous system at half-filling 
(Fig. 2(a)) and in a harmonic trap with E t = 3.28£ 
(Fig. 2(b)). For a given entropy per particle S/N the 
temperature of the trapped system is already lower than 
that of the homogeneous system at U = 0. Furthermore, 
as U is ramped up, the trapped system exhibits signifi- 
cant cooling compared to the homogeneous system. Thus 
we see that for E t = 3.28£ and any starting entropy less 
than 0.65/cb, one can obtain an AF core by adiabatic 
cooling (see Fig. 2(c)). 

We gain further insight from the profiles shown in 
Fig. 3(a,b,c). As the interaction is ramped up from 
U/t = to 8, the cloud expands and the density at the 
center decreases towards 1, characteristic of a Mott in- 
sulator (MI). This MI has a gap to charge excitations 
and thus a low entropy. On the other hand, the metallic 
state in the wings, with its low-energy spin and charge 
excitations, can act as an entropy sink. Thus, entropy is 
transferred from the Mott core to the metallic wings. 

During this process the temperature falls from T/t = 
0.53 to 0.36 ~ T/v. In the final state, the entropy density 
s(r) at the center is near the critical value for AF order- 
ing indicated by the dashed line [7]. We see the growth 
of local antiferromagnetic correlations from the nearest- 
neighbor spin-spin correlation C nn (r) = — (S r • S r+ x), 
where S r = \ J2 a pO 'a/3ci a c r p is the spin at site r. 
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FIG. 2: (a) Constant-entropy curves of a homogeneous system at half-filling [7]. There is no clear evidence for "Pomeranchuk" 
cooling as U is increased adiabatically, in marked contrast to (b). The filled symbols are QMC values for T/v from Ref. 6, and the 
dashed and dotted curves are weak- and strong-coupling asymptotic forms (see text), (b) In a harmonic trap with E t = 3.28£, 
ramping up U adiabatically produces significant cooling due to entropy redistribution. An AF state can be produced in the 
trap center even for an overall entropy per particle S/N « 0.65/cs. (c) Average entropy per particle in a harmonic trap below 
which AF order exists at the center. This is significantly higher than the critical entropy of a homogeneous system. 
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FIG. 3: Cooling by increasing interaction: Adiabatic evolution of a cloud of N — 1.3 x 10 6 particles with increasing 
interaction U for a fixed total entropy per particle S/N — 0.65/cs and trap compression Et/t — 3.28. In going from U/t = to 
8, entropy is transferred from the core to the wings (or outer shell) at r > 70, where S — J °° dr Aitr 2 s(r) remains constant. 



Our analysis shows that the adiabatic cooling in a 
trap results from entropy redistribution, and not from 
a Pomeranchuk effect in the homogeneous equation of 
state [14, 15] as discussed below. In any case, we do not 
find a significant Pomeranchuk effect (dT/dU)s < in 
DQMC, either in 3D (see Fig. 2(a)) or in 2D [16, 17]. 

Another way to cool in a trap is to use adiabatic ex- 
pansion, a standard cryogenic technique, the results for 
which are shown in Fig. 4. We see that as E t /t decreases 
from 21.93 to 3.28, the core goes from a band insulator 
to an antiferromagnetic MI. 

In Figs. 3 and 4 the open symbols used only at the low- 
est temperature (T/t = 0.36t) denote regions of the trap 
away from half- filling where the DQMC sign problem is 
significant. In this range we have used a combination of 
interpolation and results from smaller systems (for which 
the sign problem is less severe). 

We now remark on the temperature dependence of the 
double occupancy D of the homogeneous system at half- 
filling, shown in Fig. 5. As T decreases below the U ', 
D is generally suppressed due to Mott physics, so that 



(dD/dT)u > 0. At low temperature for intermediate 
[//t = 4 to 6, D shows anomalous behavior in that 
(dD/dT)u < over a range of T close to T/v- Using a 
Maxwell relation, (dD/dT) v = (dD / dS)u(dS / dT) v = 
(dT/dU) s x C/T, so that (dT/dU) s < 0, suggesting 
the possibility of "Pomeranchuk cooling" [14] by adia- 
batically increasing the interaction. However, the effect 
is smaller than predicted by DMFT and DC A [1]. This 
supports our conclusion that the "Pomeranchuk effect" in 
a homogeneous system is insignificant, as already shown 
in Fig. 2(a). 

Discussion and conclusion: 

To conclude, our most significant observation is that 
it is possible to lower the temperature of the trapped sys- 
tem by suitable adiabatic processes. Cooling results from 
entropy redistribution in a trap with the metallic wings 
acting as entropy sinks. We find that an average entropy 
per particle in the trap S/N = 0.65/c# is sufficiently low 
to produce an AF state at the center using our adiabatic 
cooling protocol. In order to go well below T/v a corre- 
spondingly lower entropy is required. 
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FIG. 5: Double occupancy D(T) for the homogeneous 
system at half-filling for U/t = 4, St = 0.125; U/t = 6, St = 
0.125; and U/t = 8, St = 0.1. The vertical red lines corre- 
spond to T N (U). 

The results for the trapped system are markedly dif- 
ferent from those for the homogeneous system. First, the 
maximum critical entropy of a homogeneous AF state oc- 
curring at U = St is 0.4&B, considerably lower than the 
average value required in a trap. Second, adiabatically 
increasing U in the homogeneous case does not lead to 
significant cooling. 

We finally discuss the implications for optical lattice 
experiments [4, 5]. Before the lattice is turned on, the ini- 
tial temperature of a trapped gas is typically Ti ~ 0.1 Tp 
where the Fermi temperature ksTp = Hujq(3N) 1 ^ 3 . For 
non-interacting fermions, an initial temperature Ti/Tp w 
0.06, within the reach of current experiments, corre- 
sponds to an average entropy per particle S/N = 0.65/cb 
in the trap. As noted above, this leads to an AF state 
at the center, which can be probed by the growth of 
nearest-neighbor spin-spin correlations. Thus, the results 



presented here imply that antiferromagnetism is achiev- 
able in optical lattices, provided that adiabaticity can be 
maintained during our cooling protocol. 

We gratefully acknowledge support from FAPERJ, 
CNPq, and INCT on Quantum Information (TCLP), 
DARPA grant no. W911NF-08-1-0338 (YLL), ARO 
W911NF-08-1-0338 (MR,NT), ARO W911NF-07-1-0576 
(RTS), and the DARPA OLE Program (RTS,NT). We 
acknowledge fruitful discussions with U. Schneider. RTS 
is grateful to Thomas Maier and Mark Jarrell for discus- 
sions regarding the comparison of DQMC with DCA. 



[1] S. Fuchs, et al, Phys. Rev. Lett, 106, 030401 (2011). 
[2] T. Esslinger, Annu. Rev. Cond. Matt. Phys., 1, 129 
(2010). 

[3] M. Kohl, et al, Phys. Rev. Lett, 94, 080403 (2005). 

[4] R. Jordens, et al, Nature, 455, 204 (2008). 

[5] U. Schneider, et al, Science, 322, 1520 (2008). 

[6] R. Staudt, M. Dzierzawa, and A. Muramatsu, Eur. Phys. 

J. B, 17, 411 (2000). 
[7] For the results in a trap we used a DQMC equation of 
state obtained from L 3 = 8 3 and St = 0.1, with a critical 
entropy sn ~ 0.32/cs at half- filling. While this slightly 
underestimates the entropy, the qualitative conclusions 
regarding entropy redistribution and cooling in a trap 
are not affected. 
[8] D. Jaksch, et al, Phys. Rev. Lett, 81, 3108 (1998). 
[9] R. Blankenbecler, D. J. Scalapino, and R. L. Sugar, Phys- 
ical Review D, 24, 2278 (1981). 
[10] S. R. White, et al., Phys. Rev. B, 40, 506 (1989). 
[11] L. Gorkov and T. Melik-Barkhudarov, Soviet Physics 

JETP, 13, 1018 (1961). 
[12] P. G. J. van Dongen, Phys. Rev. B, 50, 14016 (1994). 
[13] A. W. Sandvik, Phys. Rev. Lett, 80, 5196 (1998). 
[14] F. Werner, O. Parcollet, A. Georges, and S. R. Hassan, 

Phys. Rev. Lett, 95, 056401 (2005). 
[15] L. De Leo, et al, Phys. Rev. A, 83, 023606 (2011). 
[16] T. Paiva, R. Scalettar, M. Randeria, and N. Trivedi, 

Phys. Rev. Lett, 104, 066406 (2010). 
[17] A.-M. Dare, L. Raymond, G. Albinet, and A.-M. S. Trem- 
blay, Phys. Rev. B, 76, 064402 (2007). 



Supplemental Information for "Fermions in 3D Optical Lattices: Cooling Protocol to 

Obtain Antiferromagnetism" 



Thereza Paiva 1 , Yen Lee Loh 2 , Mohit Randeria 2 , Richard T. Scalettar 3 , and Nandini Trivedi 2 

Instituto de Fisica, Universidade Federal do Rio de Janeiro Cx.P. 68.528, 21941-972 Rio de Janeiro RJ, Brazil 
2 Department of Physics, The Ohio State University, Columbus, OH 43210, USA 
3 Department of Physics, University of California, Davis, CA 95616, USA 



PACS numbers: 



Determinantal quantum Monte Carlo: We calcu- 
late p(p,T,U/t) in the homogeneous case (no trap) us- 
ing determinantal quantum Monte Carlo(DQMC) [1, 2). 
The DQMC approach provides an exact solution on finite 
spatial lattices by retaining the full space and imaginary- 
time dependence of the fermion Green function. This is 
accomplished through the introduction of an auxiliary 
field which allows for the trace over the fermionic de- 
grees of freedom to be done analytically, the sum over 
the auxiliary field then being performed stochastically. 
At half-filling p = there is no sign problem and the 
spin correlations can be calculated down to low temper- 
atures. However, in a trap the density becomes inhomo- 
geneous and in order to capture that behavior we need 
the density away from half-filling. These calculations are 
difficult since the probability weights being used to gener- 
ate the configurations in the QMC algorithm can become 
negative. 

The QMC data is of sufficient quality that we can iden- 
tify finite-size and finite-(5r effects and perform an extrap- 
olation to L d = oo and 5r = 0. We have carried out this 
procedure for U/t = 8 and p = 0. 

The finite-^r errors in E(T) are about 0.05t (see 
Fig. 1(a)), but because they involve non-critical high- 
energy degrees of freedom, they simply scale as Sr 2 (as 
expected from Ref. 4) and can be eliminated by extrapo- 
lation (see Fig. 1(b)). The finite-size errors are smaller, of 
the order of 0.02t (see inset of Fig. 1(c)), but the discrep- 
ancy between L — 4 and L — 6 curves for 0.3 <T/t < 0.5 
is a clear signature of critical behavior, which demands 
a more sophisticated treatment. Thus, in this temper- 
ature range, we fit multiple curves E(T, L) simultane- 
ously using an appropriate scaling ansatz involving the 
known critical exponents and amplitude ratios [y « 0.70, 
a = 2 - dv = -0.11, A = « 1.52) of the 3D 

Heisenberg universality class [5], using the critical tem- 
perature from the literature [6], Tn jt ~ 0.36, and in- 
cluding a non-singular background. This allows extrapo- 
lation to the infinite-size limit, in which E(T) has a kink 
at Tn (see Fig. 1(c)). At low temperature we extrapo- 
late E(T,L) to L = oc assuming 1/L 2 scaling, and the 
result is well fit by a power law due to linearly dispersing 
antiferromagnons, E(T) = E + ciT 4 . At high temper- 
ature finite-size errors are negligible and E(T) matches 
the high- temperature series expansion from Ref. 3. 



DQMC is a complementary approach to the dynamical 
cluster approximation (DC A). DQMC gives "exact" re- 
sults for an isolated cluster (after Trotter extrapolation), 
whereas DCA simulates a cluster in a self-consistent 
bath [3]. Both require finite-size scaling. The present 
DQMC calculations have succeeded in going to low tem- 
peratures even below Tn and densities away from half- 
filling. In contrast, the present implementation of DCA 
does not include an order parameter, which precludes its 
use below Tn- Our scaling analysis suggests a critical en- 
tropy sn ~ 0.4, which is consistent with Ref. 3 to within 
uncertainties. 

Thermodynamic relations: From DQMC, one can 
calculate the density with respect to half-filling p = p— 1, 
energy per site and average double occupancy D as 
functions of /i, T, and U. From these quantities it is pos- 
sible to construct the grand potential per site ^(/i,T, U) 
in multiple ways, for example: 



Q(fjL,T) 



U 



J /2 



dp p\ 



tt(p,T) = -Tln4 + T / d/3 (E-pp). 
Jo 



(i) 



(2) 



One can thence derive formulas for all other thermody- 
namic functions, such as the entropy per site: 



dp 



dp 



dT 



(3) 



S (p,T) =0+ r 

J — o 

s(p,T)=\n4 + f3(E-pp)- d/3 (E — pp). (4) 

Jo 

Justification for LDA: The data presented here are 
obtained by simulating homogeneous (V t = 0) lattices, 
and then obtaining results for the confined system us- 
ing the local density approximation (LDA). Past com- 
parisons with DQMC simulations with Vt ^ suggest 
that the LDA provides high accuracy for the local pro- 
files of the filling, spin, compressibility, and entropy even 
in 2D [7], and hence should be even more accurate here 
in 3D. 
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FIG. 1: Elimination of 5r- and finite-size errors: (a) Energy per site E(T) for various system sizes L 3 and imaginary 
time discretization steps 5r. (b) E{5r) for L — 6 at various T, showing 5r 2 scaling. Most error bars are smaller than the 
linewidths. (c) E(T) extrapolated to 5r = and L = oo using suitable fitting forms at low, medium, and high temperatures. 
The inset shows a close-up; dashed red and blue curves are fits with scaling forms near criticality (see text). The accuracy of 
the extrapolated energy is about 0.02t. Extrapolated DCA results from Ref. 3 are shown for comparison. The extrapolated 
QMC curve (black) captures the critical fluctuations near the Neel temperature and falls below the DCA results. 
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